High Accurate Mathematical Tools to Estimate the Gravity Direction Using Two Non-Orthogonal Inclinometers

This study provides two mathematical tools to best estimate the gravity direction when using a pair of non-orthogonal inclinometers whose measurements are affected by zero-mean Gaussian errors. These tools consist of: (1) the analytical derivation of the gravity direction expectation and its covariance matrix, and (2) a continuous description of the geoid model correction as a linear combination of a set of orthogonal surfaces. The accuracy of the statistical quantities is validated by extensive Monte Carlo tests and the application in an Extended Kalman Filter (EKF) has been included. The continuous geoid description is needed as the geoid represents the true gravity direction. These tools can be implemented in any problem requiring high-precision estimates of the local gravity direction.


Introduction
The direct measurement of gravity using inclinometers has a variety of practical applications. Tilt sensors, such as accelerometers or inclinometers, are the hardware of choice for measuring a planet's gravitational field. While the gravity is perturbed on a global scale by other massive bodies (e.g., Moon, Sun), local perturbations to the gravity vector's pointing direction depend on nearby mass distributions and anomalies, or a lack thereof. A detailed review of various kinds of solid tilt sensors with comparisons between the most typical tilt sensing techniques can be found in [1].
From a sensory perspective, inclinometers are commonly used to measure inclination angle magnitudes and structural deformations. Measured values are typically provided as a percentage or angular deflection with respect to a level reference surface, whose nominal orientation is perpendicular to the local gravity vector. Other classic inclinometer use cases include measuring the slope gradient during activities such as tunneling, de-watering and excavation, as well as those that require monitoring the integrity of the ground around a structure.
In land surveying and mapping, an inclinometer can provide a rapid measurement of the slope of a geographic feature [2] or be used for analyzing the magnitude, rate, direction, depth, and type of landslides [3]. Continuous monitoring of landslides is important for early warning purposes. In general, inclinometers can be used for continuous monitoring of any movement, providing valuable information for landslide autonomous geometry reconstruction [4]. In the oil and gas industries, inclinometers are used to measure the strike and dip of geologic formations [5] in order to detect oil and mineral deposits. In forestry, tree height measurements can be made by an inclinometer using standardized methods [6]. Major artillery guns may have an associated inclinometer used to facilitate the aiming of shells over long distances [7,8]. Permanently installed inclinometers are employed at major earthworks such as dams [9,10] to monitor the long-term stability of the structure. Inclinometers can also be used in robotics [11] and to measure bridge deflections [12] for bridge safety evaluation. Specifically, the authors of [13] developed and implemented a 1.
The analytical tools to best estimate the gravity direction. These tools consist of: (a) an exact estimate of the gravity direction when inclinometer measurements are affected by zero-mean Gaussian noise, and (b) the associated gravity covariance matrix. These tools can be used in static and filtered dynamic scenarios to accurately estimate the actual gravity direction and to model the local subtle gravity variations.

2.
A continuous mathematical model to describe the true gravity direction deviation with respect to the latest (1984) revision of the approximated World Geodetic System (WGS-84) ellipsoidal model. Two deflection correction models are proposed for the north-south and east-west deflection variations. These models are expressed using combinations of a set of N linearly independent orthogonal surfaces that are derived from Chebyshev orthogonal polynomials of the first kind.
These tools and models find relevance in a variety of gravity measurement tasks and associated use, especially in filtered systems, such as the one developed for geo-technical monitoring in [22].

System Definition
The system under consideration in this study is shown in Figure 1, which describes the axes of the two non-orthogonal inclinometers, [I x , I y ], and an orthonormal inclinometer reference frame, {x, y, z}. The I x -axis of the inclinometer frame is assumed to be coincident with the x-axis reference frame, and the plane defined by the two inclinometers axes is also coincident with the [x, y]-plane. The relative orientation of the two inclinometer axes deviates from the nominal/optimal orthogonal displacement by an angle ε. This angle, which takes into account the small deviation from orthogonality due to mounting errors, can be estimated during lab calibration and, possibly, re-calibrated later to account for variations caused by vibrations and/or thermal expansions.
The gravity direction in the [x, y, z] frame is where the third component enforces the unit-norm constraint. The inclinometer measurement angles, ϑ x and ϑ y (see Figure 1), allow the estimation of the gravity direction in the reference frame.

Gravity Direction Estimation
This section contains the covariance analysis of the gravity direction estimated by two identical inclinometers, as shown in Figure 1. The analysis is performed for the generic case of a non-orthogonal mounting, where the inclinometer axes differ from orthogonality by the angle ε.
Let the two measured angles, ϑ x and ϑ y , be affected by Gaussian noise where µ x and µ y are the mean values and σ 2 is the measurement variance, identical for the two inclinometers. The mean values, µ x and µ y , constitute the two main components of the EKF state vector. The quantities ε and σ 2 are initially estimated by lab calibration methods. However, vibrations and/or thermal expansions will likely change their values over time.
In this case, these parameters should be re-calibrated and their values re-estimated by adding them as new variables to the EKF state vector (see the Discussion section). The angles, ϑ x and ϑ y , can be represented by a zero-mean Gaussian angle, δ, as To estimate the gravity direction, Equation (1) requires the evaluation of E{cos ϑ x } and E{cos ϑ y }. Dropping the axis-specific subscript, we can write, E{cos ϑ} = E{cos(µ + δ)} = cos µ · E{cos δ} − sin µ · E{sin δ}.
To compute E{cos δ}, let us expand it by Maclaurin series [23]: Reference [24] provides the following identity: Substituting Equation (4) into Equation (3), the exact expectation of the cosine of the zero-mean Gaussian angle, δ, is
This result allows us to find the expected values of the first two components of the gravity direction:ĝ g y = E{g y } = e −σ 2 /2 cos µ y cos ε − cos µ x tan ε . Figure 2 provides a geometrical interpretation of why the expectation of the cosine is lower than the cosine of the mean when the measured angle is affected by zero-mean error. If the inclinometer is not biased, then the error of the inclinometer can be described by a small cone whose axis is coincident with the inclinometer axis. This implies that inclinometer measurements affected by errors will fall within this cone. If the true angle with respect to the gravity direction is ϑ, then all measurements providing the true angle belong to the surface of a cone whose axis is the gravity direction and aperture is the angle ϑ. This cone intersects the measurement error cone by splitting it into two spherical areas that, as shown in Figure 2, are identified as "A" and "B." All measurement angles falling in the "A" area provide an angle greater than the true value, while all the measurement angles falling into the "B" area provide an angle lower than the true value. Since the "A" area is greater than the "B" area, a zero-mean inclinometer will provide an angle greater than the true value. This means that a zero-mean inclinometer provides biased angles. This phenomenon is generated by our three-dimensional spherical world, and only vanishes under two situations: (a) when the measurement angle is perfectly ϑ = 90 • , and (b) when the inclinometer standard deviation (dictating the aperture of the measurement error cone) is σ = 0 • (idealized inclination sensor). This bias effect is actually very small, even for poorly accurate commercial inclinometers. However, since the inclinometer axis direction is normally placed perpendicular to the gravity direction, this effect is even smaller, because the inclinometer works around the nominal orthogonality displacement of ϑ = 90 • .

Gravity Direction Covariance Matrix
The gravity direction covariance matrix is To compute the P(1, 1), P(1, 2), and P(2, 2) terms of this matrix, the following estimations are required: The term E{cos ϑ x cos ϑ y } can be expressed as because cos ϑ x and cos ϑ y are statistically independent, while the computation of the E{cos 2 ϑ} terms requires the estimation of E{cos 2 δ} and E{sin 2 δ}. We have The E{cos 2 δ} and E{sin 2 δ} expectations can also be computed by Maclaurin series: and, similarly, Therefore, and, specifically, Equations (9) and (10) allow us to obtain the following expectation expressions: Using Equations (6) and (11), the term, P(1, 1) = E{g 2 x } −ĝ 2 x , of the covariance matrix can be written as This expression can be simplified into the following more compact form: To derive the analytical expression of the term, P(2, 2) = E{g 2 y } −ĝ 2 y , we make use of Equations (7) and (13): After some simple manipulations and simplifications, the previous equation can be written in the more compact form Note that, as a sanity check, if the two inclinometers were perfectly orthogonal, meaning if ε = 0, then Equation (15) becomes formally identical to Equation (14).
As for the term P(1, 2) = E{g x g y } −ĝ xĝy , we have which can be simplified to This equation highlights an unexpected result: the covariance term, P(1, 2), is not a function of µ y . This result, which has been numerically verified by extensive Monte Carlo tests, remains with no explanation. Note also that, if the two inclinometers were perfectly orthogonal (ε = 0), then we would obtain P(1, 2) = 0.

Computation of the P(1,3), P(2,3) and P(3,3) Terms via Covariance Law
The analytical expressions for the three missing terms of the covariance matrix, P(1, 3), P(2, 3) and P(3, 3), were not found. However, a highly accurate estimation of these terms can be obtained using the covariance law [25,26], sometimes called the error propagation law.
The covariance law connects the two covariance matrices that are expressed in terms of two different (but equivalent) parameterizations by the Jacobian, J , of the transformation between parameterizations. In this specific case, the covariance law connects the 3 × 3 covariance matrix of the gravity unit-vector, P, with the 2 × 2 covariance matrix written in terms of the cosine of the measurements, C. This relationship consists of the identity where the analytical expressions of the covariance matrix terms , can all be computed using Equations (5), (9) and (10), while the Jacobian of the transformation between the two parameterizations is The first four elements of the Jacobian can be computed from Equation (1). Their expressions are ∂g y ∂ cos ϑ y = − 1 cos ε while the last two elements can be computed using the chain rule as Finally, the expression of the Jacobian of the transformation is where the expressions of the two partials, ∂g z ∂g x and ∂g z ∂g y , are provided by Equation (18) and Equation (19), respectively. The covariance law can actually be used to estimate all the elements of the covariance matrix. However, this law has been used here only to estimate the three missing elements of the covariance matrix, P(1, 3), P(2, 3) and P (3,3). The main reason is that the covariance law is based on a linear Taylor expansion of the transformation and is therefore an approximated method , even if it provides a highly accurate estimation of the covariance matrix. On the contrary, the terms P(1, 1), P(2, 2) and P(1, 2), provided by the Equations (14)- (16), are analytically exact.
It is easy to derive them by setting Q = C(1, 2) J(3, 1) + C(2, 2) J(3, 2). In doing so, the expressions of the last three terms of the covariance matrix can be written as

Covariance Analysis Numerical Validation
The numerical validation of the covariance matrix was performed by Monte Carlo tests using the R2016b version of MATLAB software . All of the tests were carried out using two identical inclinometers affected by zero-mean Gaussian error with standard deviation 1σ = 0.1 deg., and with axes differing from orthogonality by the angle ε = 5 deg. A set of 10 7 Monte Carlo tests were performed on 1000 gravity directions, uniformly generated in a cone with 30 deg. (π/6 rad. ) aperture around the nominal gravity direction g = {0, 0, −1} T . The uniformly distributed gravity directions were generated as where U [a, b] indicates a uniform distribution in the [a, b] range. The gravity direction components, (g x , g y , g z ), and the two angles, (ϑ x , ϑ y ), generated by this distribution have the histograms reported in Figure 3. Since the three terms P(1, 1), P(1, 2) and P(2, 2) are mathematically correct (meaning, not approximated), the Monte Carlo numerical validation tests were performed simply to quantify the accuracy of the three elements of the covariance matrix, P(1, 3), P(2, 3) and P(3, 3), whose expressions were obtained by the covariance law.
The top three plots of Figure 4 show the histograms of T(1, 3), T(2, 3) and T(3, 3), while the bottom three plots show the histograms of the differences with the three terms estimated using Equation (20). The results of these Monte Carlo tests highlight the accurate estimation of the approximated terms, P(1, 3), P(2, 3) and P(3, 3), computed via covariance law. The values of the terms T(1, 3) and T(2, 3), numerically estimated by the Monte Carlo tests, fall within the [−10 −6 , +10 −6 ] range. The differences with the terms P(1, 3), P(2, 3), estimated via the covariance law, are in the [−10 −9 , +10 −9 ] range, a difference better than three orders of magnitude with respect to the numerically estimated values. The same accuracy error level is experienced for the P(3, 3) term with respect to the numerically estimated T(3, 3) term, whose values fall within the [0, +10 −6 ] range.

Extended Kalman Filter (EKF) Validation
The variable to estimate is represented by the gravity direction, i.e., by a threecomponent vector. However, since the third component of the gravity vector, g z , is not independent (actually, it is derived from g x and g y ), then a reduced two-state EKF estimator must be developed and used. In order to quantify the accuracy gain obtained using the closed-form expression of the covariance matrix, a static scenario was selected to validate the EKF estimator. This choice does not depend on the dynamics under consideration or the accuracy of the adopted dynamical model. The EKF summary is provided as follows: • State model: Observation model: • Initialization: Given initial state vector, x − 0 , and initial covariance matrix, P − 0 , and assuming no process noise, Q k = E{w k w T k } = 0 2×2 , the measurement's covariance is The initial state, which is given by Equation (1), can be computed using the first measurements,θ x andθ y , while the state covariance matrix is where P k (1, 1), P k (1, 2) and P k (2, 1), are provided by Equation (14), Equation (16) and Equation (15), respectively. Specifically, they can be written in terms of the state vector, x − k , Predictor: Since the problem under analysis is static, • Corrector: where the Jacobian, J h (x k ), of the observation is given by where

Numerical Tests
The true direction of the gravity in the inclinometer reference frame has been selected as The inclinometer measurements are corrupted by a zero-mean Gaussian noise with a standard deviation of σ = 0.1 deg., which is a common value for commercial inclinometers. In addition, the deviation from orthogonality was selected as ε = 2 deg.  The bottom left and right plots show the absolute errors for the g x and g y components, respectively. These error plots highlight, thanks to the analytically exact expressions of the covariance matrix terms, a very fast convergence of the EKF estimator.

Continuous Gravity Description
The irregular Earth shape (and mass distribution) has presented significant challenges to surveyers and cartographers throughout history. The early map makers considered the Earth to be a perfect sphere prior to Isaac Newton positing an oblate spheroid shape model. The axial-symmetric ellipsoid description of Earth, which has been formalized as the World Geodetic System (WGS-84) [27], is reliably used as the standard shape model for most of today's GPS applications. The WGS-84 is defined by estimates of the Earth's equatorial (R E ) and polar (R P ) radii, and a flattening factor ( f ) Describing a point on the Earth's surface can be done by Cartesian coordinates (x, y, z) or by spherical coordinates, using geocentric (θ) or geodetic (φ) latitude and longitude (λ). Figure 6 shows the two definitions of latitude for an ellipsoid-shaped model. The relationship between these two latitude definitions is [28]: The greatest angular difference between geocentric and geodetic latitudes is, for the WGS-84 model, approximately 700 arcseconds, a difference that would introduce over 21 km of position error on Earth. However, even if the WGS-84 ellipsoidal model describes the gravity direction better than the sphere model, the actual gravity direction is dictated by the geoid model [29].
The geoid was first described by Gauss as the shape that the ocean's surface would take if the only acceleration forces acting on it were Earth's gravitational field and the rotational dynamics of the planet. In reality, this is still a simplification of the surface always orthogonal to the gravity direcion. The geoid is particularly useful for describing the effect that mass anomalies (such as mountain ranges or caverns) have on the local pointing direction of gravity, which can deviate from the reference ellipsoid by tens, hundreds or even thousands of arcseconds. Figure 7 shows the main parameters defining the relationship between the geoid surface and the WGS-84 ellipsoid. P is a point on the Earth's surface that is located at altitude h normal to the ellipsoid and altitude H normal to the geoid. This point on the geoid surface is located at altitude N with respect to the WGS-84 ellipsoid. The angle between the true gravity direction (normal to the geoid and provided by the inclinometers) and the fictitious gravity direction (normal to WGS-84 reference model) is the geoid correction. It is also worth pointing out that, unlike the WGS-84 axial-symmetric ellipsoid model, the geoid cannot be described by a simple equation because it is a function of the Earth's mass distribution. However, the geometry of the geoid surface can be derived with respect to the WGS-84 model using the derivatives of the Earth's gravitational potential. The Earth's gravitational potential can be described in terms of spherical harmonics: where GM = 3.986013 × 10 14 m 3 /s 2 is the Earth gravitational constant, (r, φ, λ) are the spherical coordinates at which the evaluation is made, P m n (sin φ) are the associated Legendre polynomials of degree n, order m and argument sin φ, and C nm and S nm are spherical harmonic coefficients.
Taking normalized partial derivatives of Equation (32) with respect to the geodetic latitude and longitude, we obtain a mathematical definition for the north-south (δN) and east-west (δE) deflections of the gravity direction with respect to the WGS-84 ellipsoidal model: where γ is a scalar representing the theoretical normal gravity, whose value and clear explanation are provided in Reference [30].

Continuous Description of Geoid
The analytical complexities associated with the infinite series of Equations (32) and (33) make them unable to evaluate the gravity direction without introducing truncation error. From satellite measurement data, a geoid model can be constructed in order to generate δE and δN deflection reference databases for a specified grid of points on the Earth's surface, defined by geodetic latitude (φ k ) and longitude (λ k ). To use the geoid database, interpolation between discontinuous points is needed. The accuracy of this interpolation is obtained at the expense of carrying a larger database (fine grid geoid description). In this study, a simple continuous geoid model is derived in order to use geoid corrections in specific locations and/or in dynamic situations.
The continuous geoid model is provided by expressing the δE and δN deflections as linear combinations of a set of N orthogonal surfaces: where the α k and β k are the unknown coefficients and the S k (x(φ), y(λ)) are a set of orthogonal surfaces satisfying These orthogonal surfaces are built using orthogonal polynomials, such as Chebyshev polynomials, T k (x), which are orthogonal in the x ∈ [−1, +1] range. This means that, when dealing with a specific geographical range of interest, φ ∈ [φ min , φ max ] and λ ∈ [λ min , λ max ], a mapping between the geographical coordinates [φ, λ] and the Chebyshev variables [x, y] must be defined. The most simple mapping is the linear Using Chebyshev orthogonal polynomials, the orthogonal surfaces can be defined as where any two orthogonal surfaces cannot be associated with the same [i, j] pair of Chebyshev polynomials.
Chebyshev orthogonal polynomials of the first kind can conveniently be computed in a recursive way. Starting with T 0 (x) = 1 and T 1 (x) = x, the subsequent terms are computed by the recursive relationship, The problem of fitting the geoid data grids, δE(φ i , λ j ) and δN(φ i , λ j ), consists of two over-determined linear systems: where the total number of the geoid database points is N p > N, for the problem to admit a solution. The previous equation describes two linear systems that can be written in matrix form: where the two coefficient vectors, α and β, can be computed by least squares: This continuous description of the geoid is applied in the next subsection to a geographical region where the geoid gravity direction deviations are particularly significant.

Example of Continuous Geoid Description: Himalaya Region
The results of a numerical example of this procedure are shown in Figure 8   The associated geoid database consists of a 10 × 10 grid (N p = 100 points), spanning a region of 2000 km × 2000 km. The number of orthogonal surfaces used for the continuous geoid model was N = 55. This is also the number of unknowns (size of α and β vectors). The data points of the geoid grid are shown with black markers while the continuous models are indicated by the surfaces. The L 2 norms of the least-squares residuals are shown in Figure 9. The condition number of the matrix to invert (in the least-squares process), which is provided in Figure 10 as a function of the number of orthogonal surfaces, does not exceed 200.  By selecting orthogonal polynomials (instead of monomials), the condition number of the least-squares matrix to invert is reduced and, consequently, the numerical process becomes more robust and provides a more reliable numerical estimation of the unknown coefficient vectors, α and β.

Discussion
This study introduced the mathematical tools necessary for extracting a highly accurate estimation of the local gravity (plumb-line) direction from a pair of non-orthogonal inclinometers. The statistical analyses performed herein can be used to quantify the performance of static tilt sensors, since the measured direction of gravity is known to be normal to the reference geoid surface. The next section provides an example use case of these tools for high-precision continuous gravity direction estimates under static operating conditions.
The static Stellar Positioning System [15][16][17] introduces a framework for estimating the geographical position of a vehicle on the surface of a planetary body using a pair of inclinometers along with an atomic clock and a night-sky camera to observe celestial directions (e.g., stars, visible planets). This system, which can be used in GPS-denied scenarios or in locations where the GPS signal is not available (i.e., on Moon or Mars), assumes that a reliably accurate estimate of the local gravity direction is available, which in turn would require a geoid correction. The resulting mathematical estimation problem is nonlinear and, therefore, can be solved by iterative nonlinear least squares. The mathematical procedure requires computing the Jacobian of the gravity direction and, consequently, the first derivatives of Equation (35). The first derivatives can be obtained using the recursive equation, More accurate models will require the computation of the Hessian or even higher-order terms. The subsequent derivatives needed in these cases can be computed recursively:

Future Research Directions
The authors are planning to validate the proposed gravity direction estimation method using two non-orthogonal inclinometers by: • Including the mounting angle deviation from orthogonality, ε, as an element of the EKF state vector. This would keep the system calibrated from variations of ε caused by vibrations and/or thermal expansions. • Integrating the proposed system with an Inertial Measurement Unit (IMU) to extend the use of the static "Stellar Positioning System" to dynamical scenarios. Similar to the previous point, this will extend the EKF state vector to account for the variations in the transformation matrix between the inclinometer and IMU reference frames. • Including in the EKF state estimator , for the static use "Stellar Positioning System," the transformation matrix between the inclinometer and the night camera reference frames.
As for the continuous geoid model, the authors are planning to complete the analysis on the geoid north-south and east-west deflection estimates by: • Investigating the convergence radius as a function of the geographical coordinates. • Investigating the accuracy of the convergence as a function of the geoid database grid resolution. • Including, in the current nonlinear position estimation algorithm, a divergence identification algorithm that re-initializes the estimation process using a new random guess selected within an admissible distance bound.

Conclusions
This study introduces a new approach to accurately estimate the gravity direction using a pair of non-orthogonal inclinometers affected by zero-mean Gaussian errors. The proposed approach estimates the gravity direction and its 3 × 3 covariance matrix, with the main three terms derived analytically. The accuracy of these statistical parameters is numerically validated by Monte Carlo simulation. The proposed method is then validated using an Extended Kalman Filter, which highlights a fast convergence to a very accurate solution.
Since the gravity direction is normal to the geoid surface, a continuous geoid model for the north-south and east-west corrections is provided as a linear combination of a set of orthogonal surfaces that are derived using Chebyshev polynomials of the first kind. An example application is provided for a specific region (centered in the Himalayas) to show the continuous description of the gravity direction deviations with respect to the WGS-84 ellipsoid model. This covariance analysis, along with the continuous description of the geoid, can be applied to any static or dynamic scenarios involving the measurement of gravity using a two-axis inclinometer. Funding: This research was in part supported by a NASA grant#: 80NSSC19P1369.